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ABSTRACT 

We develop a theoretical framework to construct axisymmetric magnetic equi- 
libria in stars, consisting of both poloidal and toroidal magnetic field components. 
In a stationary axisymmetric configuration, the poloidal current is a function of 
the poloidal magnetic fiux only, and thus should vanish on field lines extending 
outside of the star. Non-zero poloidal current (and the corresponding non-zero 
toroidal magnetic field) is limited to a set of toroid-shape flux surfaces fully en- 
closed inside the star. If we demand that there are no current sheets, then on the 
scparatrix delineating the regions of zero and finite toroidal magnetic field both 
the poloidal fiux function (related to the toroidal component of the magnetic 
field) and its derivative (related to the poloidal component) should match. Thus, 
for a given magnetic field in the bulk of the star, the elliptical Grad-Shafranov 
equation that describes magnetic field structure inside the toroid is an ill-posed 
problem, with both Dirichlet and Newman boundary conditions and a priori 
unknown distribution of toroidal and poloidal electric currents. We discuss a 
procedure which allows to solve this ill-posed problem by adjusting the unknown 
current functions. We illustrate the method by constructing a number of semi- 
analytical equilibria connecting to outside dipole and having various poloidal 
current distribution on the fiux surfaces closing inside the star. In particular, we 
find a poloidal current-carrying solution that leaves the shape of the flux function 
and, correspondingly, the toroidal component of the electric current the same as 
in the case of no poloidal current. The equilibria discussed in this paper may 
have arbitrary large toroidal magnetic held, and may include a set of stable equi- 
libria. The method developed here can also be applied to magnetic structure of 
differentially rotating stars, as well as to calculate velocity field in incompressible 
isolated fiuid vortex with a swirl. 
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1. Introduction 

Structure of magnetic fields in non-convective stars is a long standing issue in astro- 
physics (Prendergast f956). In fluid non-convective stars, very general arguments prove that 
both purely poloidal and purely toroidal magnetic fields are unstable (Tayler 1973; Wright 
1973; Mar key & Tayler 1973; Flowers & Ruderman 1977). As a result, magnetic fields inside 
the star should consist of a combination of both (Prendergast 1956). An example of such 
stable (on resistive times scales) equilibrium were found numerically by Braithwaite & Spruit 
(2004). No acceptable analytical solution is known. 

In this paper we discuss a semi-analytical procedure to construct magnetic structures of 
fluid stars. We are interested in developing a procedure how to construct poloidal-toroidal 
equilibria of stars, a non-trivial mathematic problem on its own. We do not discuss the 
stability of the resulting configuration. MHD simulations {e.g. Braithwaite & Nordlund 
2006) indicate that stable configurations correspond to stably stratified stars with toroidal 
magnetic field reaching inside the star the local value of the poloidal magnetic field. Thus, 
we are looking for equilibria where the toroidal and poloidal magnetic field are, generally, of 
the same order somewhere inside the star. 

Recent attempts to construct poloidal-toroidal equilibria of a spheromak-type (Broder- 
ick & Narayan 2008; Duez & Mathis 2009) suffer from a common drawback: they require 
unphysical surface currents. (Both works construct spheromak-type constant-a equilibria, 
which are stable magnetic field configurations that naturally minimize magnetic energy given 
a total helicity Woltjer (1958). Such constructions cannot avoid surface current sheets in 
principle, since poloidal current is non-zero everywhere inside the star, see discussion after 
Eq. (7)). In contrast, our approach ensures that no current sheets are formed inside the star 
or on the surface. 



2. Magnetic equilibria: the Grad-Shafranov equation 

In MHD equilibria, the Lorentz forces are balanced by gradients of pressure and gravi- 
tational forces, 

JxB = Vp + pV$, (1) 

where $ is gravitational potential. Dividing Eq. (1) by p, taking a curl, and assuming 
barotropic fluid, p — p{p) (see Appendix A for discussion of non-barotropic fluid), Eq. (1) 
reduces to 

Vx^^O (2) 
P 
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Axially symmetric magnetic fields can be written as 

B = ^Px^* + ^'^*, (3) 

where 27rP{zu, z) is the poloidal magnetic flux and /(ro, z) is the poloidal current enclosed 
by an axially symmetric loop located at cylindrical coordinates tu,z {e.g. Shafranov 1966); 
we set the speed of light to unity. Axial symmetry and the fluid approximation (so the 
pressure is derivable from a potential and thus cannot have an azimuthal component) give 
J X B|^ oc (V/ X VP) ■ = 0, implying that poloidal current is a function of P only 
/ = I{P). The force balance equation (2) gives 

Taking into account that I = I{P), this reduces to the Grad-Shafranov equation, which in 
cylindrical coordinates become {e.g. Shafranov 1966) 

A*P = -pw''F{P) - G{P) (5) 

where 

- (I) - ^- (^) 

is the Grad-Shafranov operator, G'(P) = 4JJ', / = /(P) and F = F{P) are functions of the 
flux function only. Also note, that toroidal current depends exclusively on the shape of the 
flux functions: 

A*P 

J, = -— (7) 
w 

As a crucial physical constraint, we will assume that the current cannot leave the star. 
One can then identify three regions of space, with different requirements on the functions 
I{P) and pF{P) (c. g., Monaghan 1976, Fig. 1). First, there is the outside vacuum, where 
/ = and p = 0. Second, there is a part of the stellar interior threaded by fleld lines 
connecting to the outside. Since on these fleld lines the condition I{P) = was imposed 
outside the star, it must still be true inside; non-zero density implies that the term containing 
F is non-zero inside the star. Finally, if a flux surface is completely closed within the star, 
on it / and pF can both be non-zero. 

Since we are interested in general properties of magnetic equilibria, we make a simplify- 
ing assumption of a constant density, p = constant. ^ Redeflning the function F to include 



^Equation of state and initial conditions {e.g. isentropic or not), are important for the stability of magnetic 
configurations (Braithwaite & Nordlund 2006). In this paper we are not concerned with the stability, only 
in finding stationary configurations. The derivation given below can be repeated for any given p{r). 
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the assumed constant density, the Grad-Shafranov equation becomes 

A*P = -^F{P)w^ + G{P) (8) 

The coefficient —15/2 is introduced for numerical convenience, see a comment after Eq. (9). 
Function F is assumed to be zero outside of the star (since the new definition of F inchides 
density, which is zero outside). Below we refer to function F as the pressure function (it is 
related to gradient of presser with respect to flux function Shafranov 1966) and the function 
G as the poloidal current function. 



2.1. The mathematical problem 

Stable equilibria of fluid stars must have considerable toroidal magnetic field component. 
As discussed above, the toroidal magnetic field cannot be non-zero everywhere inside the star, 
it must be limited to a smaller region of toroidal shape (in case of axial symmetry). On the 
separatrix between the regions of non-zero toroidal magnetic field and the bulk (which has 
zero toroidal magnetic field), both the poloidal and the toroidal components of the magnetic 
field must connect smoothly, without a surface current. For axisymmetric configurations, 
the toroidal components of the magnetic field is related to flux function, while the poloidal 
component of the magnetic field is related to the gradient of the flux function, see Eq. (3) 
with / = ^(-P)- In order to avoid current sheets, on the separatrix both the flux function P 
and its normal derivative 9„P should be continuous. 

Thus, we are faced with an unusual mathematical problem: we need to jind solutions of 
an elliptical equation (8) while matching on a given boundary both the flux function ( equiva- 
lently, the toroidal B-field) and the flux function normal derivative (equivalently, the poloidal 
B- field). This makes an elliptical equation (8) is formally over-constrained, as it should 
satisfy both Newman and Dirichlet boundary conditions simultaneously. In addition, the 
functions F and G are unknown a priori and should be found as a part of the solution. 
Thus, we are trying to solve an over-constrained equation with two unknown functions, 
which themselves are functions of the solution only. 

For a given magnetic field structure in the bulk (with no poloidal current), we are looking 
for magnetic field structure inside a poloidal current-carrying toroid of a given shape. In 
principle, the shape of this toroid also is not known a priori. But this does not create 
any addition mathematical complications: for any shape of the boundary, the solutions 
inside and outside of the toroid should be matched satisfying both Newman and Dirichlet 
boundary conditions simultaneously. (As a side note, a one-dimensional problem with both 
Newman and Dirichlet boundary conditions can be solved if there exist several eigenvalues 
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for the solution of either Newman or Dirichlet problem. The four corresponding boundary 
conditions can be satisfied by choosing two integration constants and two eigenvalues. This 
method does not work for the two dimensional domain, since the shapes of the flux functions 
corresponding to different eigenvalue problems are different.) In this paper we devise a 
procedure that allows to built both the functions F and G and the flux function P inside 
the toroid for given boundary conditions at the boundary of the toroid. 



3. Matching the bulk and the toroid solutions 

3.1. Solution in the bulk 

As we discussed in §2.1, we are faced with an unusual mathematic problem of continu- 
ously matching the solution of the elliptical Grad-Shafranov equation (8) with zero poloidal 
current (/ = 0) and some distribution of the toroidal current {F ^ 0) in the bulk and un- 
known non-zero poloidal current (/ ^ 0) and some toroidal current inside a toroid. As a 
basic solution in the hulk we take a well known dipolar, purely poloidal solution connecting 
to outside dipole (Ferraro 1954; Shafranov 1966) (the procedure developed below can be 
repeated for any shape of the enclosed toroid) 

B = ^ [Szuze^ + (5i?2 - - 3z')e,] (9) 

This solution corresponds to / = and F = 1. It has a set of fiux surfaces closing in- 
side a star, see Fig. (1). Under the fluid approximation, this solution is unstable to non- 
axysimmetric perturbations (Taylcr 1973; Flowers & Ruderman 1977). 

The flux surfaces are given by solution z{zu) of (9) for < P < (25/48). (Below, for 
conciseness we set R = 1, Bq = 1.) For a given meridional plane they satisfy 

9 9 5 4 P 

There is a set of nested toroids enclosing within a star and centered on the apex point 
uj — •\/5/6 and z = 0. The ellipticity of the boundary is 1/2 (the flux surfaces are not 
ellipses, though). The first closed flux surface (saparatrix) corresponds to Ps — 1/2 and 
equation for it can be written as 



(11) 
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Fig. 1. — Structure of magnetic field for constant density. Outside solution correspond to 
vacuum dipole. Inside, there is a set of nested flux surfaces near the surface and near the 
equatorial plane. Closed flux surfaces must carry poloidal current for the structure to be 
stable. 

where x = w/R — are coordinates centered on the center line of the enclosed toroids, 

corresponding to w = (2/3)^/"^, z = (since the flux surfaces are not ellipses, the symmetry 
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axis of the separatrix does not intersect the apex) . The separatrix intersects the equatorial 
plane at points vj = a/2/3 and w = 1. The maximum value of P equals Pmax = (25/48). 
This value is higher than the maximum value for vacuum dipolar field at a given radius 
Ps = 1/2; this allows a smooth sewing of the vacuum dipole potential, which is a decreasing 
function of radius, with increasing potential at small radii of the internal solution. 



3.2. The algorithm to calculate magnetic field structure inside the toroid 

The solution (9) has a set of flux surfaces closing inside a star. Toroidal magnetic 
fleld should be conflned to this region. On the separatrix between the regions of zero and 
flnite magnetic fleld, both the flux function P and its derivative should be continuous. In 
this section we discuss how the unknown functions F and G can be adjusted to satisfy the 
over-constrained boundary conditions on the separatrix d. 

We expand the unknown functions F{P) and G{P) in terms of their argument P near 
the separatrix d, corresponding to P = Pq (subscript corresponds to quantities evaluated 
on the separatrix) 

Here F^^\Po\g) and G^^\Po\g) are values of the derivatives on the separatrix. These are 
numerical coefficients to be determined. The upper bound in the sum is determined by the 
number of initial starting points and the imposed smoothness order of the solutions at the 
points of intersection. The first terms in the series (12) are known: at the separatrix F — 1, 
(7 = 0. 

At each point inside the separatrix the function P can be expanded in Taylor series 
in a following manner. On the separatrix the flux function is constant P = Po\q, while its 
derivative is normal to the separatrix and is a known function, dn Po\g. This allows us to 
express from Eq. (8) the second normal derivative at the boundary d as some function T of 
Po, dnPo, FiPolg) and G(Po\g): 

dlP\g = J^{Po\g,dnPo\a, F(Pola), G{Po\g)) (13) 

Here F{Po\g) and G{Po\g) are unknown coefficients to be fitted. By taking derivatives of 
Eq. (8) along the normal to the separatrix we can express higher order derivatives d'^'^P 
as functions of the given Pq\q and -Polg; as well as expansion coefficients F^"^^(Pola) and 
G^"^^(Po|g). The solution inside the toroid can then be represented in Taylor expansion in 
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terms of derivatives at the boundary: 



^W = E| ^"^^o|a (14) 



where s is a distance along the normal to a given point r inside the separatrix. The values of 
can then be expressed in terms of unknown coefficients F\q and 9^ G\q, similarly 

d 

to Eq. (13) for d^P\g. An explicit example of the implementation of this procedure is given 
in Appendix B. 

Next, we impose a condition that values of P extrapolated from different points on the 
boundary converge smoothly on a chosen set of points inside the toroid. Though generally 
this convergence is not guaranteed, we can consider extrapolation from a finite number of 
points on the boundary and limit the expansion to the appropriate number of derivatives 
F"^{Pq\q) and G"^{Po\q) to ensure smoothness of P at the convergence points to a given 
order. We will require smoothness of P up to second order; this ensures the absence of 
current sheets. Thus, matching of the expansions of the flux functions from different points 
gives the coefficients ^^"^^(1/2) and G^'^\l/2). This determines simultaneously both the flux 
function P and functions F and G. Since this procedure can be extended to the arbitrary 
expansion order, it determines the functions P, F and G to arbitrary precision. 

The procedure described above is implemented starting from three particular points 
on the toroid corresponding to the symmetry axes of the separatrix: at the points where 
the separatrix intersects the magnetic equator 2; = 0, r = ^y2/3Re^ and r = Re^, as 
well as at the point where magnetic fleld on the separatrix is parallel to the equator, r = 

(2/3)^/^i?en7 + -y/(5 — 2y/6)/3ez. Thus, the point where the integration curves intersect is 
zu = (2/3)^/^ and z — 0. Note, that the point where magnetic fleld on the separatrix is 
parallel to the equator has different cylindrical radius than the apex point of the enclosed 
cylinders; enclosed flux surfaces are not ellipses. This is true both for the solution (9) as 
well as the numerical solutions found here. 

Since the procedure described above is not unique, we expect various internal structures 
of the enclosed toroid. We have two unknown function, F and G: choosing one of them 
determines the other. This will create a set of various equilibria. Below we describe two 
procedures: first expansion of functions F and G in terms of P up to a given order, and 
secondly specifying F and solving for G. As is shown in §3.3, simultaneous expansion of 
functions F and G in terms of P (as opposed to arbitrary specifying one of them) selects 
a particular class of solutions, when the toroidal current and the form of the flux functions 
remain the same as in the case of zero poloidal current. (Recall that the toroidal current 
depends only on the shape of the flux function, Eq. (7).) 
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3.3. Expansion of F{P) and G{P) in terms of P - Po\g 

For three starting points on the separatrix, there are five matching conditions: equaUty 
of flux functions (two conditions), equaUties of first and second radial derivative, and zero 
derivative along the direction normal to the equatorial plane. We expand the functions F 
and G up to the third order in P — Pq\q: boundary conditions require G = and F = 1 on 
the boundary; then five expansion parameters are determined from the matching conditions 
leaving one free parameter. The free parameter normalizes the overall strength of the poloidal 
current inside the toroid. 

Results of calculations are presented in Fig. 2-3. Somewhat surprisingly, in case of non- 
zero poloidal current G and non-constant pressure function F ^ 1, the algorithm leaves 
the shape of the flux functions the same as for zero poloidal current, P{r, z) ~ Po{r, z), see 
Eq. (9) and Fig. 2. At the same time, the functions F and G on the right hand side of Eq. (8) 
are adjusted to match the zero poloidal current case, — (15/2) F{P)-uj'^ + G{P) ~ — (15/2)ro^, 
Fig. 3. 

Thus, the toroidal current and the shape of the poloidal magnetic fleld remains the 
same. In Fig. 3 the current-carrying flux functions closely resemble the current-free one, 
with precision of ~ 10^"^. We remind that this calculation is based on only three extrapolation 
points. In comparison, when the pressure function F is imposed, §3.4, the distortion of the 
form of the flux function is of the order of unity. 

3.4. Other solutions 

Next we find the fiux function P and coefficients G"^{Po) assuming that inside the toroid 
F is a given function of P. 

3.4.1. Constant F^l 

Let us find the flux function P and coefficients G"*(Po) assuming that inside the toroid 
F — ^ 1. We expand solutions to 5th order derivative near the separatrix in order to 
match expansion at the equatorial plane to the second order (three coefficients), match it to 
expansion along z axis, and condition of vanishing radial magnetic fleld at the equator. We 
investigate solutions as functions of the parameter ^. 

First, for = 1 the procedure described above closely reproduces the solution (9), see 
Figs. 4 and 5. Thus, if we keep the function F the same in the bulk and inside the toroid. 
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Fig. 2. — Structure of enclosed flux function P for the case when functions G and F are 
expanded in terms of the flux function at the separatrix. Upper row: Value of the flux function 
in the equatorial plane within the submerged toroid (upper left panel). At this scale different 
curves overlap. Upper right panel: difference between the poloidal current-carrying solutions 
and poloidal current free solution (9) in the equatorial plane. Lower row: Value of the flux 
function along the vertical symmetry axis of the submerged toroid, w = (2/3)^/^-R (low left 
panel). Lower right panel: difference between the poloidal current-carrying solutions and 
poloidal current free solution (9) along w = {2/3y^^R. 

the procedure described above converges to the solution G = 0. 

For ^ 7^ 1, there is non-zero poloidal current and, correspondingly non-zero toroidal 
magnetic field. Figs. 4. For large positive or large negative ^, the topology of the magnetic 
field inside the toroid changes, see Figs. 4, and 5 and 6. For ^ > 1.99 there forms a new set 
of flux surfaces centered close to the axis of the toroid where the poloidal field reverse. For 
—4.1 < ^ < —3.4 there is an off-centered set of flux surfaces where the poloidal field reverse. 
For ^ < —4.1 additional off-centered set of flux surfaces forms. For large absolute values of 
1^1 the sizes of the internal toroids increase, while no other change of topology occur. 

The plot of /(-P) indicates that when solutions extend to values of the flux function 
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Fig. 3. — Structure of current and pressure functions for the case when functions G and F 
are expanded in terms of the flux function at the separatrix. Left panel: Function G = All'; 
Center panel: current function I{P); Right panel: function F{P) — 1. This figure illustrates 
the value of the right hand side terms in Eq. (8) for different values of the amplitude of the 
poloidal current. 

large than, approximately .6 (in dimensionless units R = Bq = 1) or smaller that 0.4 for the 
different sign of F, the toroidal field, oc 2//r approaches poloidal field (= 1/2 at equator). 
(In case of / = the maximum value of P is 25/48.) The maximum value of P ~ 0.6 or 
minimum P ~ 0.4 is reached (approximately) for |^| > 10. Thus, we expect that for |^| > 10 
the resulting configuration is stable. 

Some of the found configurations have poloidal magnetic field reversing inside the sub- 
merged toroid. Though such configurations have not so far been seen in simulations, we 
hypothesize that this may be related to initial conditions used in simulations. 

3.4.2. F = l + a{P-l/2) 

Solutions corresponding to F = ^ 7^ 1 has a disadvantage that the toroidal current 
density experiences a jump on the separatrix. If we impose a smooth variation of toroidal 
current at the separatrix, e.g. in the form F = 1 + a{P — 1/2), the behavior of the solutions 
qualitatively remains the same. Fig. 7. 

4. Conclusion 

In this paper we develop a procedure to construct poloidal-toroidal equilibria of fluid 
non-convective stars. This involves a mathematically unusual procedure to find the distri- 
bution of unknown poloidal and toroidal currents, which depend only on the shape of the 
magnetic flux function, and the shape of the magnetic flux function itself; all these quantities 
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Fig. 4. — Shapes of the flux function P for given F = ^1 . Upper Left Panel: analytical 
solution Eq. (9). Upper Right Panel: numerical solution for ^ = 1. In this case the numerical 
solution closely matches the analytical solution Eq. (9); this vindicates the fitting procedure. 
Lower Left Panel: ^ = 10. Lower Right Panel: ^ = —10 

must match the shape of the separatrix and specified Dirichlet and Newman boundary con- 
ditions on it. The procedure we propose uses expansion of the current functions in powers 
of the flux function to a specified order, which depends on the number of points used to ex- 
trapolate the flux function. Generally, this procedure is not unique; for example, given some 
pressure function F it will determine the poloidal current function G, or vise versa. This is 
qualitatively consistent with results of numerical simulations Braithwaite & Nordlund (2006) 
which show that depending on the initial conditions a variety of final states are achieved. 

One particular configuration stands out among the all possibilities. It corresponds to 
simultaneous expansion of two functions F and G in terms of the value of the flux function 
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Fig. 5. — Current as function of flux function P, P in the equatorial plane, P along the line 
w = (2/3)^/*^ as function of z , Bz in the equatorial plane for different values of parameter 
^ = 1, ±10 (except in the I{P) plot where ^ = 1, ±5, ±10. The middle line in the plots (b-d) 
is both the analytical solution (9) and numerical solution for ^ = 1 (the difference between 
them is tiny on this scale). 

P on the separatrix. In this particular case, the shape of the poloidal current-carrying flux 
surfaces remains the same as in the case of no poloidal current: newly found functions F 
and G add up to produce the toroidal current of the poloidal current-free configuration. 

An important feature of our solutions is that there no current sheets, neither on the 
surface of the enclosed toroid nor on the surface of the star. On the other hand, we do not 
test for the stability of the resulting magnetic structures. This could be done, in principle, 
by minimizing magnetic energy at a given helicity for magnetic field structure inside the 
toroid. Results of numerical simulations (Braithwaite & Nordlund 2006) indicate that stable 
magnetic field configurations require both stable stratification and typically have toroidal 
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magnetic field of the order of the poloidal magnetic field. Our procedure, in principle, allows 
arbitrary ratios of toroidal and poloidal magnetic fields, and we expect that some of the 
found magnetic field configurations to be stable. 

One of the main limitations of our approach is that we assume barotropic equation of 
state. In contrast, in stably stratified stellar interiors buoyancy forces play an important role, 
dominating typically over magnetic forces. Still, even in non-barotropic fluid the poloidal 
current is still function of P (Appendix A), so that the toroidal magnetic field is still limited 
to a set of fully submerged fiux surfaces. On the separatrix of the regions with vanishing and 
finite toroidal magnetic field, again, both the fiux function and its derivative should match. 
We leave this problem to future considerations. 

The method developed here can also be applied to rotating stars when the magnetic 
axis and the rotation axis are aligned. In this case the Grad-Shafranov equation has a form 
similar to (8) with F — > F + zu'^QQ', where ^{P) is the angular velocity of rotation, which 
is constant on a given flux surface (Chandrasekhar 1956). The new function 0.{P) can be 
expanded near the saparatrix in a same way as the pressure F and the current function G. 

There is a hydrodynamical analogue to the solutions presented here, describing an iso- 
lated fluid vortex with a swirl. Stationary magnetic configurations are related to velocity 
field of a time-independent flow of an incompressible fluid, with the substitution v — > B 
(Shafranov 1966). In case of a steady, axially symmetric fluid motion, the velocity potential 
satisfles the same equation as (5), the so-called Bragg-Hawthorne and/or Squire- Long equa- 
tion (see, e.g. Lamb 1975, Eq. 165.13). Our method then describes the swirling velocity of 
a core of an overall non-swirling vortex (see also Moffatt 1969). 

I would like to thank Jonathan Braithwaithe, Martin Kruczenski, Daniel Phillips, An- 
dreas Reisenegger, Hendrik Spruit and Dmitri Uzdensky. 
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A. Non-barotropic EoS 

In non-barotropic fluid Eq. (2) becomes 

„ J X B Vp X Vp , , , , 

V X = „ Al 

P P 

Axial symmetry and the fluid approximation still imply / = /(P). The force balance equation 
(Al) gives 

VP . v^Zlili: = VpxVp 

writing p = p(s, p), where s is entropy and assuming 

A*P + 4/ J' 

— J— = P(.,P), (A3) 
w^p 



This preprint was prepared with the AAS L^T[t;X macros v5.2. 
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we find 

VPxVF^^-fvPxVs=(^l] (A4) 



ds \ds J ^ p 



If we assume tliat distribution of density is splierically symmetric, p = p{r), we find 

Tfiis form of F should be used in the Grad-Shafranov equation (8), which makes it an 
integro-differential equation for P {cf. Villata & Ferrari 1994, Eq. 20). Thus, in case of 
non-barotropic fluid the poloidal current is still function of P only, I — I{P), plus there will 
be an extra term in the Grad-Shafranov equation, Eq. (A5). Thus, one still encounters the 
same problem, that on the boundary of the enclosed toroid with toroidal magnetic field both 
flux function and its derivative should match those in the bulk. 



B. Explicit example of expansion of P near the point w — 1, z — 

Let us illustrate the procedure of finding the flux function inside the toroid starting from 
a particular point on the separatrix, zu = 1, z = 0. In this case the normal coincides with 
the radial direction. We know the value of Po\q = 1/2 and its derivative 9„ Pq\q = —1/2 
at this point. In addition the flux function is zero on the separatrix, G{l/2) = 0, while the 
pressure function is unity, F{l/2) ~ 1. From (8) we find d^P : 



d^P 1 



dlP = -dlP + ^ - -^^'F{P) + G{P) (Bl) 
w Z 



On the boundary it is equal to 

1 15 



dlP\a -—2- ^^"nm - G(l/2) 



13 



(B2) 



Taking a derivative of Eq. (Bl) with respect to vo we find 



dlP\, = -18 + f F'(l/2) + ^G'(l/2) (B3) 
Thus, expansion of the flux function near the point vo — 1, z — Q reads 
P 4 _ _ 1, _ ^ , (_18 + H^'(V2) . 1g'(V2)) <^ . ... (B4) 

It can be continued to arbitrary high order. Matching of the expansions of the flux functions 
from different points then gives the coefficients F("*)(l/2) and G'('")(l/2), thus determining 
simultaneously both the flux function P and functions F and G. 
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Fig. 6. — Qualitative structure of the poloidal flux surfaces closing inside a star for F = ^ 7^ 1. 
For —3.4 < ^ < 1.99 {Panel (a)) the direction of the poloidal magnetic field inside the toroid 
is of one sign. For very large, ^ ^ 1 {Panel (b)) there is a reversal of magnetic field close to 
the axis of the toroid. For large negative, ^ ^ —4.1 {Panels (c,d)) the reversal of magnetic 
field occurs inside two separatricies not aligned with the main axis. Thick lines indicate 
separatricies, where magnetic field reverse direction. These figures are illustrations only. 



0.85 0.90 055 lOO 0.85 050 0.95 1,00 0,00 0.05 0.10 0.15 

r/R r/R z/R 

Fig. 7. — Form of the flux functions and magnetic field for a given F = 1 + a{P — 1/2). 
Qualitatively, solutions remain the same as in the case F = ^ ^ 1. Left panel: values of P at 
the equatorial plane. Center panel: values of magnetic field at the equatorial plane. Right 
panel: values of P at zu = (2/3)^/^. Qualitatively, the structure of the solutions remain the 
same as in Figs. 5-6 



